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ABSTRACT 

Whereas the higher-order versions of the finite element method (the p- and 
hp-versions) are fairly well established as highly efficient methods for monitoring 
and controlling the discretization error in linear problems, little has been done to 
exploit their benefits in elasto-plastic structural analysis. In this paper, we discuss 
which aspects of incremental elasto-plastic finite element analysis are particular 
amenable to improvements by the p-version. These theoretical considerations 
are supported by several numerical experiments. First, we study an example for 
which an analytical solution is available. It is demonstrated that the p-version 
performs very well even in cycles of elasto-plastic loading and unloading, not 
only as compared to the traditional h-version but also in respect to the exact 
solution. Finally, an example of considerable practical importance - the analysis 
of a cold-worked lug - is presented which demonstrates how the modelling tools 
offered by higher-order finite element techniques cam contribute to am improved 
approximation of practical problems. 
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1. Introduction 


Elasto-plastic analysis is of increasing importance in all fields of engineering 
sciences. However, real-life elasto-plastic problems are difficult to solve because they are 
non-linear in nature. Approximate numerical methods are required for solving practical 
problems of elasto-plasticity. The last two decades have witnessed a great deal of research 
for finding approximate solutions to elasto-plastic problems by the finite element method 
(FEM). For an excellent review of the techniques of elasto-plastic modeling in the 
framework of the finite element method, we refer to [1]. 

Almost all research work on elasto-plastic finite element analysis was based on the 
traditional h-version of the finite element method. The theoretical basis of the h-version is 
explained in [2]; see, e.g., [3] for details on the application of the h-version to elasto- 
plasticity. Whereas a lot of attention has been focussed on realistic material models, the 
equally important issue of the accuracy of numerical modeling, including questions of 
discretization errors and convergence of the error in energy norm, remain largely 
unaddressed (see, however, [4], [5]). 

During the last decade, techniques for monitoring the discretization error in finite 
element approximations have been studied thoroughly, and new strategies for an efficient 
minimization of the error in energy norm have been developed. In the p-version of the 
FEM, the polynomial degree of the finite element approximation is increased to improve 
the accuracy of the solution, rather than following the traditional approach of mesh 
refinement (h-version). The hp-version combines mesh refinement and increasing 
polynomial degrees. 

The theoretical basis and convergence properties of the h-, p-, and hp- versions of the 
FEM for linear elliptic problems have been well established. For a survey of the state-of- 
the-art of the p- and hp-versions of the finite element method based on the displacement 
formulation, the reader is referred to [6]. By all measures of performance, the higher- 
order methods provide superior accuracy (with respect to the discretization error) than 
their h-version counterpart for a large class of linear problems including linear elasticity. 

The performance of the p-version for elasto-plastic stress analysis problems has not 
been investigated until very recently. To the authors' knowledge, the first numerical 
experiments which include incremental elasto-plastic material models in the frame of 
higher-order finite element methods have been reported in [7]. However, at that time no 
further discussion or evaluation of the efficiency of the p-version for this class of 
problems has been given. The first author of the present study has developed an 
experimental finite element code named FEASIBLE (cf. [7]) which includes p-version 
capabilities and provides elasto-plastic material laws as well as multi-stage analysis. This 
implementation has been employed for the numerical examples included in sections 5 
and 6. 
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Recently, the applicability of the p-version for solving elasto-plastic problems, using 
the deformation theory of plasticity, has been investigated numerically in [8]. This study 
shows that, even though the deformation theory is, strictly speaking, restricted to cases 
where the principal stresses remain proportional to each other throughout the plastic 
process, both theories yield very similar results in many cases. However, the problem 
categories amenable to the deformation theory do not include unloading and cyclic elasto- 
plastic processes. 

Therefore, the present paper focusses on the incremental theory of plasticity (for an 
introduction, see, e.g., [1], [9]). The paper is organized as follows: We identify those sub- 
tasks of the elasto-plastic finite element procedure which determine the computational 
performance and accuracy of the numerical analysis. Subsequently, we discuss in which 
steps of the elasto-plastic analysis we would expect improvements when switching to the 
p-version. It is shown that many of the advantages of the p-version in purely linear 
analysis carry over to elasto-plastic problems. A numerical study of the loading and 
unloading of a thick-walled cylinder under internal pressure is presented. For this 
problem, an exact solution is available. The analytical results are then compared with the 
finite element solutions obtained by the p-version and the h- version. We demonstrate that 
the p-version method performs very well, also as compared with the h-version. 

Finally, some unique features of the p-version finite element code are illustrated 
which permit realistic modeling of practical engineering problems. The finite element 
solution to a lug which undergoes cold-work expansion is provided and compared with 
the approximate closed-form solution often used in practice by engineers. This example 
demonstrates that many cold-working problems demand a finite element analysis when 
other than very crude accuracy requirements are imposed. Furthermore, it shows that the 
p-version of the FEM makes the analysis of such problems very convenient and easy as 
well as reliable and robust. 



3 


2. Finite element approximation 

Once a specific material model has been selected and the geometrical and mechanical 
idealization of the structure has been determined, the problem of structural analysis is cast 
in the form of a purely mathematical problem; The finite element method is a tool for 
solving this well-defined mathematical problem approximately. The primary aim of the 
finite element analysis is to find as good an approximation to the exact solution of the 
mathematical problem as can be obtained with a given amount of engineering work and 
computer resources. The issue of how well the selected mathematical model corresponds 
to the physical reality is to be treated completely independently of the issue of the 
numerical accuracy of the finite element solution. In the present paper, we focus 
exclusively on the errors introduced by the finite element discretization, not on the 
physical modeling issue. 

The finite element method based on the displacement formulation starts by 
minimizing the functional of the potential energy 

n(u) = jB(u,u)-F(u) (1) 

over a finite-dimensional subset 5 of the space S(Q.) of functions which have finite 
strain energy over the domain Q and fulfill the displacement boundary conditions. Here, 
the bi-linear form ^ B(u, u) denotes the strain energy corresponding to some displacement 

field u and F(u) the work of the applied loads. Thus, an approximation u FE e S(Q), 

minimizing n(fi) over S(Q), is obtained, whereas the exact solution u ex minimizes n(fi) 

over S(Q) . 

Consequently, for an assessment of the efficiency of a specific finite element scheme, 
the magnitude of the discretization error in the energy norm is monitored as a function of 
the number of degrees of freedom N used in the finite element approximation u FE . The 
error in energy norm is defined by 

l|e|| E = Jn(u ex - u FE ) . (2) 

In the present study, the p-version of the FEM is applied in the form of a uniform p- 
extension, increasing the polynomial degree of the approximation uniformly on the entire 
mesh. As far as linear analysis is concerned, the p-version is distinctly superior to the h- 
version for problems which have an exact solution that is analytic throughout the domain. 

In such cases, exponential convergence rates (i.e., ||e|| £ =ce~^) can be achieved 
asymptotically (N — » with simple finite element meshes, whereas the best 

convergence rate of the h- version is only algebraic (i.e., ||e|| £ = cN^). 
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Even for problems which include a finite number of points where singularities occur 
in the stresses, the p-version converges (algebraically) at least twice as fast as the h- 
version. Resorting to the hp-extension, the exponential convergence rate can be recovered 
even for these singular problems. Convergence rates significantly better than algebraic 
can be obtained in the pre-asymptotic range (i.e., for practical N) by performing the p- 
extension on a geometrically graded mesh in these cases, cf. [6]. 

In elasto-plastic analysis, n(u) is a non-linear function of u so that the properties of 
the FEM for linear analysis do not carry over to elasto-plasticity directly. However, in the 
incremental theory of plasticity, we solve a sequence of linearized problems with the 
same geometry, but strain-dependent material properties and loads. The final solution is 
obtained as the sum of those of the quasi-linear steps. Therefore, the overall performance 
of the non-linear analysis will be determined by the smoothness of the exact solutions of 
the stepwise problems. 

The smoothness properties of the exact solution of purely elastic problems are well 
known and can generally be determined a priori from the geometric shape and the 
boundary conditions; on the other hand, little has been done yet to inquire how much the 
introduction of elasto-plastic material behavior changes the smoothness properties of the 
corresponding purely elastic solution. 

However, when deciding which finite element method is more appropriate when 
solving an elasto-plastic problem, the question of smoothness is of little relevance 
because the higher-order methods perform significantly better for problems in any of the 
two categories of smoothness discussed above. This important fact is often overlooked. 

Furthermore, it is worth noting that the elastic-plastic material is still characterized by 
continuous strains so that elasto-plasticity does not introduce any kind of an abrupt 
"material interface" inside the elements. Such a situation might cause the p-version to 
yield oscillatory answers. However, it does not occur in small-strain elasto-plasticity since 
there is always a continuous change in the material properties, even though the one- 
dimensional ideal elasto-plastic stress-strain relationship intuitively suggests some kind 
of "discontinuity". Even upon unloading and reloading, the strains and consequently the 
material properties (directly depending on the strains) remain continuous. 

All contained plastic flow problems fall into the category of completely continuous 
strains and, consequently, the corresponding displacements can be expected to be at least 
not much more unsmooth than those resulting from purely elastic problems. 

As far as uncontained plastic flow is concerned, the underlying mathematical problem 
changes from the elliptic to the hyperbolic type, including possible bifurcations, 
multiplicity of the solution, and the like; such problems are strictly beyond the limits of 
small-strain elasto-plastic analysis for both the h- and p-versions of the FEM. They 
cannot be solved by methods based on the principle of minimum potential energy. Of 



5 


SHSsrSaS “ r • - 

mathematical model has toT^d £ to w* * diffcrcnt 

approximation. ° d> ^ mm W* o/ ./Writ trfemimr 

ppilSs^ts 

prolonged incremental computation with many load steps. 

daealv“om tU the MtlT^T nUm “ CaUy dle pointwisc 9 -aIity of strains computed 
n , , y h f eIement approximation of the displacement field for linear 

poblems. This is exactly the technique to will be employed an <Ll l^Z 

HOI itotifi 1 , ‘° eteimme 46 current material Properties. Hie numerical investigation 
[10] identifies various patterns of convergence and shows that in all cases and for a^ide 

rhtloTkT™ore ti0n , faC, ° rS H <a T 0111 "' 55 atai “> 46 

me h version. Furthermore, it is evident to veiy small errors in the energy norm are 
usuaily required to achieve sufficient accuracy of the pointwise strains, even for quke 

smooth solutions. Such global accuracy levels are impractical for an h-version approach 
due to the slow convergence of this method. approach 

. ^ ™ nsiderations indicate that the p-version of the FEM will be 

beneficial for elasto-plastic computations, and there is no indication that the p-version 

mig give rise to any new numerical problems not encountered in the h-vers^on. Our 
numerical examples strongly support this. 



6 


3. Numerical aspects of J 2 plasticity 

In the present paper, we are dealing with the application of the p-version of the FEM 
to problems from metal plasticity. The computations are based on the ideal elastic-plastic 
J2 flow theory, using the von Mises yield law. Further assumptions include small 
displacements and small strains. The incremental formulation of the von Mises yield law 
starts with the following assumptions: 

The total strain increment can be decomposed into a purely elastic and a purely plastic 
part: 


dz = dz el +(k pl 

We assume the engineering definition of the strain vector. 

The yield criterion is independent of the hydrostatic pressure and of the third invariant 
of the deviatoric stress tensor: 

F = 3J 2 -G 2 y . (4) 


J2 denotes the second invariant of the deviatoric stress tensor and o y the uniaxial yield 

stress of the material. No strain hardening is assumed. Only stress states such that F < 0 
are permissible. Once the stress path reaches F = 0 in a point, plastic strains will develop 
in that point upon further loading. 

During a plastic loading increment, the stress state is confined to the yield surface 
given by eq. (4): 

dF= 


dF 

9a 


dc r = 0. 


(5) 


In equation (5), dc denotes the stress increment vector and the superscript t indicates 
transpose. 


The direction of the plastic flow is given by (normality rule) 

dF 1 


dz pl = dk 


9a 


( 6 ) 


Here, dk is a proportionality factor and a denotes the stresses in vectorial notation. 
Using Hooke's law for the elastic part of the strains, 

ck el = D~ l da , 


(7) 
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where D is the linear-elastic material matrix, formulas (3), (5), and (6) can be combined 
to determine the elasto-plastic stress increment 


D 


da = D ep dz = i 


D — 


'<^1^1 V 

3a 3a 


' 3 F 

3 a 


-if 


D 


3 r 

3a j 


cfe. 


( 8 ) 


corresponding to a total strain increment cfe. Formula (8) defines the elasto-plastic 

material matrix D ep . This relationship is needed in a finite element program to determine 

the tangential stiffness matrix for a Newton-Raphson method or a similar scheme. 
Furthermore, it can be used to integrate the constitutive law on the integration point level 
of the finite element method in case an explicit algorithm is to be used for that. Whereas 
the correctness of the tangent stiffness matrix has only an influence on the efficiency of 
the non-linear solver, the integration of the constitutive law has a direct impact on the 
accuracy of the final results. 

Much attention has been dedicated to the accurate integration of the constitutive law 
during the last two decades (cf. [11], [12], [13], [14]). A comprehensive scheme of error 
control and quality assurance in the FEM has to include this aspect, which is, however, 
independent of the finite element discretization technique. With a view to practical 
applications, and especially with respect to the complicated plastic laws which become 
increasingly common in rock mechancis and concrete modeling, the authors of the 
present study favor a fully general implicit scheme similar to the one presented in [14], 
However, at present, only the traditional tangent-radial return method has been 
implemented. This simple scheme can be considered sufficiently accurate for the von 
Mises yield criterion without hardening provided that small steps are used for the 
integration of the constitutive law. 

However, the most challenging computational issue specifically related to the use of 
J2 flow rules is created by the condition that the plastic strains correspond to an 
incompressible mode of deformation (note that it follows form equations (4) and (6) that 
the plastic strain increments have no volumetric component): Starting from a 
compressible elastic behavior, the material becomes progressively more incompressible 
during the elasto-plastic deformation process as the amount of plastic strain increases as 
compared to the elastic strain. This will eventually lead to a severe locking problem in the 
traditional h-version of the finite element method if the displacement formulation is used. 
The problem of incompressibility locking in J 2 plasticity has already been studied in [15], 
and various schemes have been devised to avoid the locking, either by "reduced 
integration" or by introducing an independent approximation of the hydrostatic pressure. 
Today, most approaches use the mixed finite element method resulting from the latter 
idea, entailing all the disadvantages and additional complications inherent to mixed 
methods. 
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On the other hand, the problem of locking due to incompressibility has been studied 
extensively for the p-version of the finite element method (see [16], [17], [18]), and it has 
been shown that no locking effects occur when the polynomial degree p is greater than 4. 
No special precautions or mixed methods have to be introduced to obtain accurate results. 

Our numerical experiments will demonstrate clearly that locking due to 
incompressibility does not occur in the p-version analysis of the elastoplastic problem. 

4. Implementation 

In many details, the implementation of an elasto-plastic material law in a p-version 
FEM code is identical to a standard h-version implementation. Therefore, we summarize 
only briefly the main features of the implementation used for the present study: 

First, we carry out a purely linear elastic analysis and perform a uniform p-extension, 
starting at p=l (piecewise linear approximation). The maximum polynomial degree 
available is p=8. The convergence of the strain energy is monitored. Based on the 
assumption that the error in energy norm converges algebraically, a very reliable estimate 
of the exact energy can be obtained by extrapolation from three different p-levels (cf. [6]). 
Throughout the present study, the space of hierarchical tensor product trial functions has 
been used (product space). 

To avoid errors in the linear solution by approximate geometrical mapping, we always 
use the blending-function method (see [6]) to describe curved boundaries accurately. 
Therefore, we can use very large elements even in the presence of curved boundaries. 

When the error estimate for the purely linear solution reaches a user-defined 
threshold, the non-linear computation starts. There is no guarantee that the non-linear 
solution will be of the same order of accuracy as the linear solution, but we rather expect 
the nonlinear solution to be less accurate. We have already seen that, e.g., 
incompressibility effects which are not originally present in the system will be introduced 
by the plastic material law, and those effects will slow down the convergence of the error 
in energy norm. Therefore, it seems reasonable to impose very restrictive accuracy 
requirements on the linear solution. In our practical computations, as a rule of thumb, we 
set the desired accuracy of the linear solution to about 10% of what we would consider 
reasonable for a purely elastic analysis. 

The next step is to determine at which fraction of the total imposed load the system 
will start to yield. The remaining load is divided into a number of load increments, and in 
each of these, the Newton-Raphson method with a consistent tangent predictor is applied. 
The Newton-Raphson iterations are stopped when both the relative magnitude of residual 
force vector and of the change in the incremental displacements in either the Euclidean or 
the maximum norm are smaller than a prescribed threshold. In each load increment, a 
check for elastic unloading is performed in every integration point. 
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Finally, stress results are available only in the integration points since the incremental 
elasto-plastic law has been evaluated only there. To obtain stress results elsewhere, a 
least-squares approximation for each of the components of the stress vector is computed 
tor each element individually. For this least-squares fit, we employ the trial functions 
corresponding to the polynomial level used in the finite element approximation of the 
displacements. 

5. The thick-walled tube under internal pressure 

First, we study an example for which an analytical solution is available, the thick- 
walled tube under internal pressure. The analytical solution for both the loading and the 
unloading - including elasto-plastic re-loading in the reverse direction - has been given in 
[19]. It is one of the few elasto-plastic two-dimensional problems for which an analytical 
solution even for the unloading is availabe, and therefore we use it as a basis for the 
evaluation of the numerical quality of the h- and p-versions of the finite element method. 

However, this problem is also of considerable practical interest since the pressurized 
tube may be used as a model for the cold-working of holes, e.g., in aircraft engineering. 
Rich et al. [20] suggest to use the analytical solution for the analysis of a wide range of 
cold-worked holes, and in our practical application example, we will study how well this 
solution is suited for a cold- worked lug. 

Cold-working is a process which is used to produce favorable stresses that 
significantly reduce the effects of stress concentration around a hole. In this method, an 
oversized and tapered mandrel, sometimes with a lubricated sleeve, is pulled through the 
hole to be cold-worked. Upon the removal of the mandrel, the hole is surrounded by a 
region of residual compressive stresses . This method is generally used for holes at the 
time of aircraft manufacture, as well as during service, to increase the fatigue life of 
structural components. 

5.1 Analytical solution 

The analytical solution for the problem of the tube under internal pressure has been 
given in [19] and applied to the cold- working process by [20]. The following assumptions 
are adopted in the analysis: 

• elastic-perfectly plastic material model 

• plane strain situation 

• von Mises yield condition 

• incompressibility of both the elastic and plastic zones (no analytical exact solution 
exists for compressible elastic deformation). 

In the application to cold-working, the mandrel is assumed to remain purely elastic. 
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Notation. Consider a tube with an internal radius a and an external radius b. Young's 
modulus Ep, Poisson ratio v p and uni-axial yield stress a y . Denote the modulus of 

elasticity for the mandrel by Efr, its Poisson ratio by v 6 and its radius by a+I. The 
difference / between the radii of the mandrel and the hole is called interference. 

Upon inserting the mandrel into the tube, a plastic zone of radius p is created. This 
radius is defined by the following implicit equation (if no solution exists for the equation, 
the interference is simply too small to create a plastic zone): 



0.52 


W+i-g' 

a b j 


♦is*#. 

a E p 


= 0, 


and the pressure between the mandrel and the hole can be shown to be 



H +1 ~¥ 


(9) 


( 10 ) 


2cr, 


If -j=-ln—, then the entire cylinder becomes plastic, and this case is of no interest 

Vo CL 

for us. 


The circumferential stresses around the hole are given by: 


Oe = i 


__y_ 
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y 

V 3 


21 n — + 1+-21 

P 

P* P 2 


b 2 + r 2 


a< r< p 
p < r < b. 


( 11 ) 


Upon removal of the mandrel, the tube undergoes an elastic unloading. If 




a 


-\±— Zy*)’ no reverse yielding occurs, and the residual circumferential stresses 
are given by: 
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a < r < p 
p < r < b. 


( 12 ) 


Otherwise, compressive yielding occurs in the zone a < r < p R , where p R is given by the 
implicit equation: 


01 a , p R PV 3 . 

2 In l + + — — = 0. 


Pr 


2a„ 


(13) 


In this case, the circumferential stress is: 
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1 
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a<r<p R 
p fi < r < p 
p < r < b. 


(14) 


This problem is suitable for both the incremental and deformation theories of plasticity 
because it is axisymmetric and the stresses are proportional. 

Note that eq. (6) in [20] is not correct which renders the stresses in the compressive 
plastic zone around the hole after unloading. The correct solution is found in [19]. 

5.2 Numerical Examples 

The problem described in the previous section is solved by the p-version as well as 
the h-version of the FEM. The finite element solutions to this problem will demonstrate 
the capabilities of the p-version and its accuracy. 

Comparison with the exact solution 

We study a tube with an internal radius a=1.5in and an outer radius b=4.0in subjected 
to a cold- working process. A plane strain situation is assumed. A pressure of P=62,000psi 
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is applied to both the inner side of the tube and the outer face of the mandrel. We 
determine the corresponding interference level as the sum of the radial displacements of 
both bodies. Thus, an interference 7=0.03399in is obtained. The material properties of the 
tube are defined by £ p =10,000,000psi, the uniaxial yield stress a y =58,OOOpsi, and 
Vp=0.49. We select this value of the Poisson ratio in order to be close to the exact 
analytical solution which assumes elastically incompressible material. The displacement- 
based finite element schemes cannot be applied if v=0.5 exactly. The mandrel's material 
properties are taken to be ££=30,000,000psi, and v^=0.3. All finite element computations 
are based on an elastic-perfectly plastic von Mises yield law. 

Exploiting symmetry, only one quarter of the domain was analyzed. The radial 
displacements are constant along the perimeter. For each of the two versions of the FEM, 
two different meshes have been analyzed. The purely elastic mandrel has been studied 
only in one simple p-version mesh, imposing extremely high accuracy requirements. 


D C 




Figure 1: Meshes used in the p-version analysis of the tube under internal pressure. 

In Fig. 1, the meshes used for the cylinder in the p-version are shown. We use a 
uniform mesh of two elements and a graded mesh of six elements for the representation 

of the tube. The boundary conditions are u y = 0 on edge A-B, u x = 0 on edge C-D, and 

t n = —P on A-D, where t n denotes the normal traction. The other boundary conditions are 

traction-free. The p-level is increased on each element until the estimated global error in 
energy norm is decreasing below 0.1%. 

For the h-version, we use meshes consisting of 70 8-noded elements and 140 8-noded 
elements, respectively. Both meshes were graded in the radial direction such that the 
biggest element was ten times bigger than the smallest. A mesh consisting of 300 8-noded 
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elements was also analyzed to make sure that the mesh of 140 elements was fine enough 
in the circumferential direction. The three different h- version meshes are shown in Fig. 2. 
The results obtained with the 300-element mesh are practically identical to those of the 
140-element mesh, so only the latter mesh is used in our examples. 



70 elements 140 elements 



300 elem. 

Figure 2: Meshes used in the h-version analysis of the tube under internal pressure 

Furthermore, we also analyzed the mesh with 70 elements, employing 9-noded 
elements instead of the 8-noded elements. The trial function space associated with the 9- 
noded elements corresponds to the product space at p=2 in the p-version. 

The h-version computations were carried out with the commercial finite element 
solver ADINA [21]. The AD ENA system offers a special pressure-displacement element 
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for incompressible materials, based on the mixed finite element formulation [22], This 
element was used in the analysis. 



Figure 3: Circumferential stress in pressurized tube, loading case. Results of p-version 
FEM analysis and exact solution. 

In Fig. 3, the exact analytical solution of the circumferential stress (<Tq) for the 
loading case is compared with the p-version solution. The h-version solution with 8- 
noded elements and the exact solution are displayed in Fig. 4. Note the very good 
agreement between the exact and the finite element solutions in the whole range, except 
at the interface between the elastic and plastic zones. At this location, the h-version 
solutions with the 8-noded elements tend to oscillate, and the quality of the results is 
questionable, whereas the p-version solution is much smoother and follows the exact 
solution more closely. The p-version merely produces smoothed stress results close to the 
elastic-plastic interface. It is clearly visible that the elasto-plastic radius (p) can easily be 
determined from the p-version solutions, but not so well from the solutions obtained by 
the h-version with 8-noded elements. 

However, the h-version solution improves considerably when we employ the 9-noded 
elements instead. These elements have an additional "internal mode" trial function. The 
corresponding results are displayed in Fig. 5. No oscillations arc present any more, even 
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in the coarser mesh. The computer time for this analysis is somewhat greater than for the 
8-noded elements. It is interesting that, even in the h-version, obviously some degree of 
"p-extension" and switching to a different space of trial functions is required to capture 
the solution of this problem properly. 



Figure 4: Circumferential stress in pressurized tube, loading case. Results of h-version 
FEM analysis with 8-noded elements and exact solution. 
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Figure 5: Circ umf erential stress in pressurized tube, loading case. Results of h-version 
FEM analysis with 9-noded elements. 

Upon removing the internal pressure, the tube first undergoes an elastic unloading; in 
the last stage of the unloading, reverse yielding occurs in the vicinity of the inner surface 
of the tube. The overall performance of the h-version as well as the p-version, compared 
to the exact analytic solution, is presented in Fig. 6 and 7, respectively. 



Circumferential Stress (psi) 
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Figure 6: Circumferential stress in pressurized tube. Residual stress after cold-working. 
Results of p-version FEM analysis and exact solution. 
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Figure 7: Circumferential stress in pressurized tube. Residual stress after cold- working. 
Results of h-version FEM analysis and exact solution. 

In Figs. 8, 9 and 10, we present "zoomed" plots of the circumferential stresses in the 
reverse yielding zone. The same phenomenon as in the loading case is observed at the 
new elasto-plastic interface radius p K : Both h-version meshes with 8-noded elements 
exhibit oscillations of the stress in the reverse yielding zone. Furthermore, it turns out that 
the 2-element p-version mesh is obviously unable to capture the very narrow zone of 
reverse yielding properly. Clearly, just one element in the radial direction is not sufficient. 
However, a very good solution is obtained with the 6-element mesh and the p-version. 
Also, the solution of the h-version with 9-noded elements is quite accurate. 
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Figure 8: Circ umf erential stress in pressurized tube. Residual stress after cold- working. 
Results of p-version FEM analysis and exact solution. Zoom of zone of compressive 
plastification. 

In order to verify the smoothness of the material properties in this elasto-plastic 
computation, the condition number (ratio of the largest and smallest eigenvalue) of the 
tangen tial material matrix (stress-strain-matrix) has been checked in all integration points 
during the loading and unloading process. The condition number changes very little 
during the loading and unloading and depends predominantly on the elastic material 
properties. This indicates that the influence of elasto-plasticity on the smoothness of the 
stress-strain law is very small. 




Circumferential Stress (psi) 
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Figure 9: Circumferential stress in pressurized tube. Residual stress after cold-working. 
Results of h-version FEM analysis with 8-noded elements and exact solution. Zoom of 
zone of compressive plastification. 
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Table 1: Computation times and number of unknowns for the eiasto-plastic loading 
and unloading of a tube under internal pressure (v = 0. 49). The FEA solutions were 
computed on a Silicon Graphics SGI IRIS Indigo workstation (R3000, SPEC mark=26). 
The computation times reported below do not include post-processing operations. The 
number of degrees of freedom does not include constrained degrees of freedom. 



ADINA 

FEASIBLE 

Coarse mesh 

82 s 

88 s* 


450 d.o.f. 

288 d.o.f. 

Fine mesh 

162 s 

90s** 


870 d.o.f. 

456 d.o.f. 


* m aximum p-level is 8, product space 
** maximum p-level is 5, product space 


The realistic case ofv p =Q.3 in the elastic zone 

Usually, the Poisson ratio of engineering materials in the elastic region is close to 0.3. 
In the following, we use the finite element solution to show that the analytical solution 
based on the incompressibility assumption in the elastic region may be considered as a 
good approximation for the solution with v p = 0.3. 

For this case, we select a smaller interference than before, 7=0.02in. This will 
emphasize the influence of the compressibility in the elastic zone since the extent of the 
plastic zone will be smaller than before. Employing the 6-element mesh with the p- 
version, and the 140-element mesh with the h- version, the corresponding pressure was 
calculated to be P=52,200psi. The finite element solutions for v p = 0.3 are compared to 
the analytic solution (assuming v^= 0.5) in Fig. 11. It is evident that the unrealistic 
assumption of incompressibility in the elastic region does not significantly affect the 
validity of the analytical solution for the practical case. The eiasto-plastic radius, 
however, is smaller than the one predicted by the FEM solutions. 
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Figure 11: Circ umf erential stress in pressurized tube, v = 0.3. Residual stress after cold- 
working. Results of p-version FEM analysis and analytical solution for v = 0.5. 

6. Practical application: The cold-working of an attachment lug 

Usually, as a first approximation, engineers use the closed form solution for a cold- 
worked cylinder given in [20] in real life problems, even though the geometry of the 
realistic problem is different from the cylinder. A cold-worked fastener hole located near 
a plate edge, or a lug, is usually idealized as a cold-worked tube having an inner radius 
equal to the hole radius, and an outer radius equal to the distance of the fastener hole from 
the edge. We show in the following that this approximation has considerable deficiencies 
and should be used only for fairly crude accuracy requirements. 
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Figure 12: Mesh used in the p-version analysis of an attachment lug. Asterisks indicate 
the maximum extent of the plastic zone during cold-working. 

In the case of the lug, which is no longer axisymmetric, the question of adequate finite 
modeling arises. The plastic zone is no longer bounded by a circle, and also the 
deformations of both the mandrel and the hole deviate from the purely circular shape 
Therefore, a finite element analysis imposing uniform pressure or constant radia^ 
displacements on the hole does not reflect the true situation properly. Neither the final 
shape of the deformed mandrel and hole nor the pressure distribution is known a pnon. 
Furthermore, even for the axisymmetric case of the tube (with corresponding 
dimensions), an analysis carried out with prescribed displacement equal to the 
interference (i.e. assuming a rigid mandrel) reveals considerable differences in the 
residual stresses (roughly 10% deviation) in the elastic zone of the tube when comparing 
with the deformable mandrel case. Also, the plastic radius is over-estimated by the model 
with the rigid mandrel, although the stresses in the plastic region are very close to those 
obtained with an elastic mandrel and the corresponding pressure inside the tube differs by 

only 1%. 

In the following, we show how to pursue a more realistic analysis. 

The lug with a hole of radius a=1.5in, outer radius Z>=4.0in and total length of /=12in 
was modelled by 13 elements (see Fig. 12). An interference level 7=0.02m was assumed 
between the mandrel and the lug. For the lug, we use material constants 
E- 1 0,000, OOOpsi, (?y=58,OOOpsi and v=0.3. Symmetry conditions were used along edges 
A-B and C-D to reduce the computation time. 

We assume that the influence of the elastic stiffness of the mandrel can be represented 
by a constant-stiffness distributed elastic spring interface element attached to the edge o 
the hole. The corresponding approximate normal spring stiffness is readily obtained from 
the axisymmetric linear analysis of the mandrel subject to constant pressure The 
tangential stiffness of the spring is zero (no friction). Displacements corresponding to the 
original interference between the lug and the mandrel are now prescribed on the interface 
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elements as an initial displacement rather than directly on the hole. Thus, we obtain a 
deformed shape of the hole (and the mandrel) which deviates from the circle, and a non- 
uniform pressure distribution reflecting the actual stiffness relation between the lug and 
mandrel. 

Note that the lug is constrained exclusively by the spring interface elements and the 
symmetry constraints. We consider the cold-working process in a co-ordinate system 
attached to the mandrel. The analysis simulates a cold-working procedure in which the 
mandrel is not capable of carrying any resultant force in the plane of the lug, but is 
allowed to move freely. Therefore, we obtain an equilibrated pressure distribution along 
the hole. The lug as a whole moves, but the left end of the model (attached to the airplane 
and thereby constrained) remains almost a straight line provided it is far enough from the 

hole. 

After loading, removal of the mandrel is simulated by deleting the interface elements. 
At this stage, the model of the lug has only rigid body and symmetry constraints. Traction 
boundary conditions are given by the pressure formerly carried by the spring interface 
elements. We are now solving a problem which has only traction boundary conditions 
(except the symmetry constraint), and therefore we require the tractions to be in 

equilibrium* 

Note that the two stages of the computation are characterized by different boundary 
conditions and, consequently, a different number of unknowns, and that adding or 
deleting elements in an existing model may change the smoothness of the exact solution. 
This has to be kept in mind when checking the accuracy of the corresponding purely 
linear solutions to ensure the use of reasonable models in the two-stage non-linear 

analysis. 

Four interface elements were used along arc B-C. The material properties of the 
mandrel are the same as in all previous examples, giving a normal spring stiffness 
C =38,500,000psi. Fig. 10 displays the maximum extent of the plastic zone during the 
loading as obtained by the p- version analysis with p=5 (product space), corresponding to 
an estimated error in energy norm of less than 0.1% in the elastic solutions. The 
maximum deviations of the internal pressure during cold-working from the constant value 
P=52,200psi (corresponding to the tube) were +8.7% and -2.5%. 
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Figure 13: Cold-worked attachment lug. "Circumferential" ’stresses after -unloading for 
vjSous locations. Comparison with analytical solution of cold-worked tu . 

In Fig 13 the "circumferential" stress in the lug after complete unloactog is 
compared with the solution obtained by application of the analyncal so ^° n ^ ^ ^ 0 
Hie unloading is purely elastic. Fig. 13 displays the stress at 9=00, 9=90°, and 9- ISO 

close to the outer boundary of the lug. 

This realistic problem clearly demonstrates that the simple thick-walled tube analogy 
and unloading conditions. 
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Conclusions 

We conclude that the p-version of the finite element method performs very well in 
small-strain incremental elasto-plastic analysis of structures. The p-version makes it a lot 
easier to control the discretization error of the finite element method. Error monitoring 
and control is achieved automatically, with minimum user interne non, and without 
substantial loss of computational efficiency. In some cases, we fin at e mo e s ase 
on the p-version yield clearly superior accuracy m comparison with their h-versio 

counterparts. 

As an example, the paper demonstrates the application of the p-version of the FEM to 
the cold-working problem. It is seen that the method is very well suited for this practically 
relevant problem and that reliable and realistic results can be obtained with a small 

amount of modeling work. 

In summary, the application of the p-version in elasto-plastic structural analysis can 
be recommended. 
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